Critical Dynamics of a Vortex Loop Model for the Superconducting Transition 
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We calculate analytically the dynamic critical exponent zmc measured in Monte Carlo simulations 
for a vortex loop model of the superconducting transition, and account for the simulation results. In 
the weak screening limit, where magnetic fluctuations are neglected, the dynamic exponent is found 
to be zmc = 3/2. In the perfect screening limit, zmc = 5/2. We relate Zmc to the actual value of 
z observable in experiments and find that z ~ 2, consistent with some experimental results. 
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PACS Numbers: 05.70.Jk, 74.40. +k, 75.40. Gb, 75.40.Mg 



The discovery of the short coherence length cuprate su- 
perconductors has allowed heretofore inaccessible fluctu- 
ation effects in superconductors to be probed. Beginning 
with the penetration depth measurements of Kamal et 
al. Q , and including measurements of magnetic suscep- 
tibility 0|§, resistivity HQ and specific heat H, static 
and dynamic fluctuation effects have been convincingly 
observed and accurately quantified. These measurements 
are consistent with the theory of a strongly type-II su- 
perconductor, with a weak coupling of the order parame- 
ter to the electromagnetic field, described by the 3D XY 
model coupled to a gauge field Q. 

The dynamic critical exponent, z, characterizes the 
relaxation to equilibrium of fluctuations in the critical 
regime of systems exhibiting a second order phase tran- 
sition 0J^]. In particular it relates the time scale of re- 
laxation, t, to a relevant length scale, x: t ~ x z . For 
infinite systems x is the correlation length, £. Near the 
critical point, the correlation length diverges and the re- 
laxation time tends to infinity, a phenomenon known as 
critical slowing down. In finite size scaling studies, x is 
identified as the system size L. 

The dynamic critical exponent, obtained from the 
measurement of longitudinal dc-resistivity for YBCO is 
z = 1.5 ±0.1 in finite but small magnetic fields ||. Simi- 
lar results were reported for the zero-field DC conductiv- 
ity pf]| , pl| . Frequency dependent microwave conductivity 
experiments yield z ~ 2.3 — 3.0 12 1. On reanalysis it was 
found that the data were consistent with z ~ 2 provided 
one neglected the region close to T c p3). Moloni et al. 
obtained z = 1.25 ± 0.05 at low magnetic fields H], but 
a later, more complicated analysis by these authors gave 
z = 2.3 ± 0.2. More recently, DC conductivity measure- 
ments on single crystal BSCCO samples were interpreted 
to give evidence for z ~ 2 |fL5| . In summary, experiments 
do not yet yield a consistent picture of the critical dy- 
namics. 

If the dynamic exponent were indeed z ~ 1.5 , then this 
would be surprising. Precisely this value is obtained for 
the superfluid transition in He 4 where the combination of 
second sound (a propagating mode, therefore z = 1) and 



order parameter dynamics (diffusive, therefore z ~ 2) 
lead to z = 3/2 (model E dynamics) @. In YBCO, 
however, the combination of a momentum sink arising 
from the lattice, and the Coulomb interaction destroying 
the longitudinal current fluctuations should lead to pure 
order parameter dynamics and a prediction that z ~ 2 
(model A dynamics). It is of course possible that some 
other mechanism can yield z ~ 1.5. 

To shed light on these issues the critical dynamics was 
investigated numerically by performing a Monte Carlo 
calculation of z for the 3-dimensional XY model, in the 
vortex representation (the so-called Villain model p6|), 
with and without magnetic screening fl7|| . The spin wave 
degrees were replaced by discrete vortex variables and the 
dynamics imposed was dissipative. The dynamic expo- 
nent estimated through a scaling analysis of the resistiv- 
ity calculated within linear response will be denoted by 
zmc- Surprisingly enough the exponent was found to be 
zmc ~ 1-5 when the interaction was unscreened while 
zmc ~ 2.7 in the presence of screening. Not only does 
the value zmc ~ 1-5 agree with previous results obtained 
by performing a similar analysis on the London Lattice 
Model (LLM) |ll| but also with the value of z reported in 
some of the experiments cited above. The observations 
in the computer simulations are surprising because there 
are no collective modes in the Villain model so that the 
dynamics would be expected to be purely diffusive, with 
zmc ~ 2. Nevertheless, and contrary to expectation, 
here too the system seems to support model E dynam- 
ics. Other extensive simulation studies report values of 
zmc ~ 1-5 and zmc ~ 2 depending upon the boundary 
conditions 19 1. 

The purpose of this Letter is to calculate analytically 
the dynamic exponent for the Villain model. The equa- 
tion of motion, corresponding to the Monte Carlo steps 
implemented in the numerical computation, is derived 
and analyzed near equilibrium. A scaling analysis is used 
to extract Zmc- We are able to explain the simulation re- 
sults in both strong and weak screening limits. We show 
also that the simulation results cannot be interpreted as 
providing evidence in support of the z ~ 1.5 result found 
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in some experiments, because they do not measure the 
true dynamic critical exponent: zmc z - We show how 
to relate zmc and z, and find that the result zmc ~ 1-5 
is in fact an artifact of taking the thermodynamic limit 
and the range of vortex interactions to infinity limit in 
the wrong order. The correct physical prediction from 
the simulation is z ~ 2 for any finite range of interaction, 
consistent with some observations. 

The Villain model:- Consider the XY model with a 
fluctuating vector potential a represented as lattice gauge 
theory link variables a™ = a,j — a,- : 
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where J is the coupling constant, Ao is the screening 
length, <pi is the phase of the condensate on site i of a sim- 
ple cubic lattice of size N = L? with periodic boundary 
conditions. The first sum is taken over nearest neigh- 
bors, while the second is over plaquettes of the lattice. 
The lattice spacing has been set to unity. The fluctuating 
gauge potential atj satisfies the constraint that at each 
site i, the discrete divergence vanishes: [V • a]i = 0. The 
phase degrees of freedom can be replaced by vortices by 
introducing the periodic Villain function to replace the 
cosines. Standard manipulations [^0| lead to the dual 
Hamiltonian: 



ijGij [Aq 



(2) 



where the n^s are vortex variables that reside on the links 
of the dual lattice and C?y is the screened lattice Green's 
function, 



Gij [Ao] — J 



(2^) 2 



exp[zfc • (f - fj)] 



L3 Y 2 E™[i-cos(fe m )] + V 



(3) 



The two limits that are considered in the simulations 
are the long range case, Ao — ► oo, and the short range 
case, Ao — > 0. Actually the simulations were performed 
by setting Ao = and Ao = oo in (3). The distinction 
between the limit and the actual simulations will turn out 
to be significant. In both cases the local constraint of no 
monopoles, [V • n]i = 0, is imposed. Each Monte Carlo 
move consists of trying to create a closed vortex loop 
around a plaquette. The trial state is accepted or rejected 
according to the heat bath algorithm with probability 
1/[1 + exp((3AE)] where AE is the change in energy and 
(3 = 1/ksT, with ks being Boltzmann's constant. Each 
time a vortex loop is formed it generates a voltage pulse, 
AQ = ±1, perpendicular to its plane, the sign depending 
on the orientation. This voltage fluctuation gives rise to 
an electrical resistance, R, which can be analyzed within 
linear response theory. A point that will be important 
to note here is that R depends on the average change in 
the total number of loops pointing in a given direction at 



each time step. The unit of time is normalized so that, on 
average, an attempt has been made to create or destroy 
one loop per plaquette. 

Dipole gas description:- It is known that near T c the 
static properties are dominated by the proliferation of 
vortex loops of unit strength, i.e., it is energetically unfa- 
vorable to create vortex loops of greater strength at each 
plaquette. The interaction between these vortex loops is 
spherically symmetric and so is the state in thermal equi- 
librium. As it stands, the computations above have been 
performed on what is known as the low temperature Vil- 
lain model and the critical point is obtained by looking at 
the intersection of the low and high temperature Villain 
models (for details see ref. pOfl). The physics described 
here is that of an interacting gas of dipoles, d. In the long 
range case they interact via the standard Coulomb term 
which falls off as 1/r 3 ; note that these dipoles interact 
antifcrromagnetically, and are not current loops, which 
interact via the standard ferromagnetic interaction. 

For our analysis we shall consider a cubic lattice, L 3 , 
on whose vertices reside the loop variables, U. In terms 
of the vortex variables n*j = V x as can be seen by 
writing out the components. The three components are 
each either ±1 or 0, corresponding to a clockwise, an- 
ticlockwise or absence of a vortex loop along the three 
principle directions, x, y or z. The corresponding proba- 
bilities on site i at time step s are given by -Pj"[l], P^[~ 1] 
and .Pj"[0], where a is a coordinate label. The quantity 
computed in the simulations is the total number of loops, 
N" pointing along a given direction a at time step s. 



^+i = E( p 5+i[i]-^+i[-i]) 



(4) 



To study the behavior of N", we follow the standard 
procedure of writing out the master equation for the time 
development of the probabilities and evaluating (4) [5l| . 
As previously indicated, the equilibrium state is spheri- 
cally symmetric. That is, on average, AEf s , the change 
in energy on adding a unit loop on site i at time s, is zero. 
This implies that transition probabilities for creating and 
annihilating a vortex loop are equal. The heat bath al- 
gorithm ensures that the conditions of detailed balance 
are satisfied. Furthermore, at T c , the restriction to unit 
loops per plaquette results in P?[0] = PQ[i] = P£,[— 1] in 
equilibrium. Since we are interested in small deviations 
from equilibrium, we impose a uniform perturbation, 5l a 
per site and see how it relaxes back to equilibrium. This 
implies SN a ~ L 3 8l a . To leading order the equation of 
motion reads 



d8N a 
dt 



sr 



(5) 



where the subscript denotes equilibrium, and af is 
the transition probability in equilibrium for creating the 
dipole loops. 

Scaling analysis:- Equation (5) is the basis for the scal- 
ing analysis that follows. The only relevant length scales 
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are the system size, L, and the correlation length, £. af is 
an equilibrium microscopic transition probability which 
remains finite at the critical point while (3 J^i A-B" is di- 
mensionless and scales as (L/£) 3 away from T c for finite 
systems. This follows because by definition, thermody- 
namic additivity occurs on a scale beyond the correlation 
length. While the free energy is extensive for all temper- 
atures, at T c , £ ~ L and /3£\A_E? is independent of 
L. Thus the characteristic time scale of relaxation of the 
perturbation, r, scales as 



(6) 



where [I] is the scaling dimension of the field I. 

For the long range case the binding energy is given by, 



0H 



— dj ~ ' f ^ d i ' til) (7) 



where di — fili, fi is the dipole strength of a unit loop 



around a plaquette, 



|, where = n 



and 



fij is the unit vector along ry. If U were dimensionless, 
then the energy of the system would not be extensive. To 
evaluate the dimension of I note that L 6 [l] 2 /£ 3 ~ (L/£) 3 
as required by the extensivity of the free energy. Thus 
[I] ~ L~ 3 / 2 and r ~ £ 3 L~ 3 / 2 . The dynamic exponent 
at T c , where £ = L, in this case is zmc — 3/2, which is 
consistent with the computer simulation results. 

For the short range case the binding energy is given by 



/3^(V x di) • (V x di) (8) 



Requiring extensivity, i.e. [l] 2 ^ 2 L 3 ~ (L/£) 3 , yields 
[I] ~ ^~ 1 / 2 . From (6) we get r ~ £ 5 / 2 which at T c scales 
as L 5 / 2 . The dynamic exponent is zmc — 5/2, which is 
consistent with the computer simulation results |T^ ]. 

Critical dynamics of the dipole gas model:- We will now 
derive the governing stochastic partial differential equa- 
tion that describes the long wavelength critical dynamics 
of the superconductor. Our strategy will be to first derive 
the continuum limit of the Hamiltonian (2), and then im- 
pose relaxational dynamics. We will find that the results 
for z are not the same as our results for zmc- This is 
because the Monte Carlo time step does not correspond 
to the physical time step. This is explained below. Let 
us first look at the continuum limit of the short range 
case. Reintroducing the coupling constants and the lat- 
tice spacing, a, we write the Hamiltonian Hy for the 
vortex variables as 



Ar 



H v = J I 2tt— x k) ■ (V x k) 



(9) 



Converting the sum to an integral, 



H v = -(J/a 3 ) ( 27r ^) J d r{V X 1(f)) ■ (V x 1(f)) 

(10) 

In the limit a — > 0, Ja~ 3 — » J and Xo/a — > Ao. Re- 
laxational dynamics is governed by the time-dependent 
Ginzburg-Landau equation (TDGL), which in this case 
is 



g=rj(2^A ) 2 (V 2 f 



V(V-Z)) + ?f 



(11) 



where r\ is a white noise, satisfying the fluctuation dis- 
sipation theorem with (r] a (f)) = and (rf 1 (f )rj^ (f)) — 
2Tk B TS a0 8(f - f). The TDGL equation is similar to 
the diffusion equation and is expected to yield a dynamic 
exponent of z = 2, in mean field theory, with small cor- 
rections due to fluctuations. The linearity of the TDGL 
in this case reflects the fact that only unit vortices are 
considered in the analysis. 

In the long range case, taking the continuum limit, we 
obtain 

Bv.jv.rl** (fxf(r ;»- ( r rw) (12) 

J \r—r\ 

x exp[-|f- fl/Xo] 

where the infinite self energy has been subtracted, and 
the screening length Ao is taken to be finite. To relate 
this to the dipole-dipole interaction used in our analytic 
model of the simulations, consider a cubic lattice, as be- 
fore, on whose vertices sit variables di, and take Ao = oo. 
Replacing 2n\/jl(f) — '^2 i di5(f — ft) one can perform 
the integrals over f and r', to recover the expression in 
(7). The actual TDGL equation for the long range case 
reads 



dl(f) 
dt 



FJ(2tt) 2 J 



df 



V 2 l(f) - V(V • 1(f)) 

\f-f] 

x exp [-(\f-f\/Xo)] +rj 



(13) 



Let us first take the case L — * oo with Ao finite but 
large. The dynamic exponent in this case is 2, because 
the kernel effectively renormalises the time scale in a way 
that is independent of system size. If we took the two 
limits L — > oo and Ao — ► oo in the opposite order, as 
was done in the computer simulations, the exponential 
factor would not be present, and the dynamics would be 
independent of L. Hence the dynamic exponent would 
then be z = 0. 

Nature of the long range case:- The rather curious re- 
sult of z = is obtained for the situation where the 
screening length is sent to infinity before taking the ther- 
modynamic limit. Whether the interaction is considered 
short or long range depends on with what it is compared. 
Physically the short range case describes the situation 
where Aq is much smaller than the inter-vortex spacing 
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A„ . This is indeed captured in the simulations by set- 
ting Ao = 0. Physically the long range case describes the 
situation where L 3> Ao X v . This is not captured by 
setting Ao = oo with L finite. 

Reconciliation with the lattice model simulations:- The 
critical dynamics of the lattice simulations and the con- 
tinuum analysis above do not apparently agree. We now 
will show that this is because the time step in the simula- 
tion does not correspond to the physical time step. The 
reason is that from the definition of the loop variable h, 
the net electric field at time t is E a (t) = £\ lfdP[l?]/dt. 
In the simulations, and in (4), this has been replaced 
by E Q (t A ,/c) = J2i dP[l?]/dtMC, where t M c denotes the 
Monte Carlo time and If = ±1,0 only. However in the 
long range case and in the short range case at T c , where 
£ = L, [I] depends on L. Hence the physical time is re- 
lated to the Monte Carlo time by t = tMc[l] so that the 
relaxation time is actually 

r~£ 3 [Z] 2 K. (14) 

The dynamic exponents for the lattice model are then 
z = 2 for the short range case and z = for the long 
range case, in agreement with the analytic calculation 
based on the continuum limit equations of motion. We 
see that the simulation result zmc = 3/2 or equivalently 
its corrected form z = arise from taking the thermo- 
dynamic limit and the long-range of interaction limits in 
the incorrect order. With this correction to the results of 
the simulation, the results no longer are consistent with 
those experiments reporting z ~ 1.5. 

Experimental ramifications:- In experiments performed 
on bulk superconductors one would expect the short 
range limit of the model above to apply, provided that 
the interaction range is shorter than the system size. In 
such systems, as long as diffusive dynamics for the vortex 
degrees of freedom is applicable, a dynamic exponent of 
2 is predicted by the model above. 

What then could be the origin of the behaviour z ~ 1.5, 
if confirmed, in some experiments? There are at least two 
possible avenues for further investigation into the true na- 
ture of the critical dynamics in these systems. The first 
is to seek experimental evidence for the existence of hy- 
drodynamic modes which might account for the observed 
model E dynamics in transport properties. The 41 meV 
peak observed in neutron scattering data is one possible 
candidate p2j although it does not seem to occur near 
the origin, while certain interpretations of the peak-dip- 
hump structure seen in ARPES are also suggestive of the 
existence of a collective mode in the system ^3|. A sec- 
ond possibility is to study the crossover from model E to 
model A dynamics as the effective coupling of the con- 
densate with the electromagnetic field tends towards zero 
(equivalently one can study the crossover by sending the 
plasmon gap to zero). 

In conclusion, we have explained the simulation results 
for the critical dynamics of the superconducting transi- 
tion in zero field, and shown that in fact they are consis- 
tent with expectations based on the TDGL. An extension 



of this analysis to two dimensions will be presented else- 
where. 
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